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ABSTRACT 

Pair annihilation of heavy stable particles that occurs in the early universe is 
reconsidered including the off-shell effect not properly taken into account by the 
conventional Boltzmann equation approach. Our new calculation of the time evolu- 
tion shows that the off-shell effect prolongs the freeze-out, always with a larger final 
relic abundance. The final yield (number density/temperature 3 ) is insensitive to the 
effective coupling for the annihilation and of order 10~ 8 x (M/l GeV) 1 ^ 3 , with M 
the heavy particle mass, if the coupling is not too small. 
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The dark matter problem in cosmology makes it worthwhile to accurately cal- 
culate the relic abundance of weakly interacting massive particle (WIMP) in modern 
particle theories. In most past works the annihilation rate has been computed using 
a thermally averaged Boltzmann equation 0J. But estimate based on the on-shell 
Boltzmann equation is dubious at very low temperatures, and importance of the 
off-shell effect for the unstable particle decay has been discussed in [fjj. We extend 
these works to the annihilation process such that it can be applied to estimate the 
relic abundance of the cold dark matter. 

Let us first state our main result and its consequences. The yield Y, the number 
density per comoving volume, is defined relative to (cosmic temperature) 3 ; Y = 
n/T 3 . Introducing the inverse temperature x = M/T for the time (M being the 
mass of the heavy particle), the yield follows the evolution equation, 

dY 



— = - r] (x) (r 2 -r cq (x) 

(2) 

The temperature dependent parameter r)(x) is roughly the thermally averaged (on- 
shell) annihilation rate (cr a v) n divided by the Hubble rate at the temperature T, 
and is, at T = M, of order 10~ 3 A 2 m p i/M with A a small dimensionless coupling 



and m-pi the Planck mass. Here d = y47r 3 iV/45 , related to the effective number of 
particle species N contributing to the cosmic energy density. 

The equilibrium abundance Y eq (x) is usually given by the thermal number density 
of the ideal gas form and becomes at low temperatures Y eq (x) ~ (x/2tt) 3 ^ 2 e~ x for 
particles without any conserved quantum number. This is very small due to the 
Boltzmann factor e~ M l T at x — M/T 1. One may well question whether this 
formula is still valid at the low temperature relevant for the freeze-out. We shall 
show below that the off-shell effect modifies this to 

The first new term represents a power dependence on the temperature (ex T n+1 ), 
with 5 containing a coupling factor, for instance O[10~ 6 ]x (4-body coupling) 2 , in a 
simple model, the boson-pair annihilation model we later explain. The power n = 
for S-wave annihilation and n — 1 for P-wave annihilation. In accord with the higher 
angular momentum n, the temperature dependence of the rate is r](x) = i]/x n+2 . 
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We plot in Fig.l the solution to the new time evolution equation, ([[]). We took the 
parameters relevant to the boson-pair annihilation model, eq. (|I~9"D . For comparison 
a result based on the Boltzmann equation without the 6 term is also shown. 

The analytic estimate of the freeze-out temperature [BJ, important for the relic 
abundance, is as follows. One may consider with a good precision that the yield 
follows the equilibrium abundance Y e<i prior to the freeze-out temperature. For the 
simplicity of this estimate we assume that the freeze-out occurs when the new 5/x n+1 
term dominates over the exponential one in Y eq (x). The freeze-out temperature Tf 
is then determined by dY eq /dxj = —r]Y^ q /x 7 j +2 , since after this epoch the inverse 
process is almost frozon, along with dY/dx = — r]Y 2 /x n+2 . Integration of this 
equation gives the final yield Yq as T — > 0. When Yf is very small and xj 1 <C 1, 
Y faYf, as is usually the case. 

In terms of the model parameters, 

One sees that the freeze-out yield Yf is insensitive to the coupling, since in simple 
models the coupling factor cancells in the ratio 5/rj, depending on the mass ratio 
as 5/r)(xM/m p \. In the boson model Y s « 2.4 x 10~ 8 iV 1/6 (M/l GeV) 1/3 . The 
contribution to the mass density in the present universe is po = MYfT$ with T ~ 
3 K. Numerically, 

po « 4.1 x 10 4 N 1/6 (M/l GeV f 13 eV cm- 3 , (5) 



using the ratio 5/r] for the boson-pair annihilation model. Note for comparison that 
the on-shell Boltzmann equation gives Tf/M ~ 1/lnr/, Yf ~ (lni]) n+2 /ri . 

The closure mass density of order 1 x 10 4 eV cm~ 3 thus excludes the WIMP mass 
range far above 1 GeV. In Fig. 2 we show the contour map of the present mass 
density in the parameter (5/rj,yS) space, assuming the mass relation between 5/r] 
and M/iripi valid for the boson model. The parameter region in which the new term 
dominates over the Boltzmann suppressed term in Y eq (x) is 0.73 d 1 ^ 3 > (In r/) 16 r\~ 2 l 3 . 

In some parameter region of supersymmetric theories A in the boson model is 
replaced by ~ g 2 (M/M) 2 . Here g is a generic gauge coupling constant, and M 
is a generic mass of exchanged superparticles for the four-body process, usually 
much larger than the LSP (lightest supersymmetric particle) mass M. With g 2 = 
Ana w 0.1, a favored case of M/M = 0.1 and M = O[10 GeV] is within the off-shell 
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dominant region given by M < 10 2 GeV , and is excluded by the overclosure from 
Fig. 2. It is thus wise not to neglect the off-shell efffect in calculation of the LSP 
abundance. Due to the P-wave mixture and detailed numerical factor not known at 
present, a realistic calculation in SUSY is called for. 

We then turn to a sketch of how the equilibrium abundance Y eq (x) above is 
derived. A complete derivation involves a fair amount of technical details, and in 
this short note we only give the essence of derivation, relegating technical details 
and proofs of the statement to our longer companion paper H| . 

For simplicity we take an interaction hamiltonian density of the form, |<^ 2 x 2 , 
where ip is a heavy bosonic field and \ a lighter bosonic field. The dimensionless 
coupling A must be less than unity for our result to be a good approximation. The 
annihilation channel ipif — > XX is related to the scattering channel (fx ~^ fX by 
a crossing symmetry. It is important to recall that a finite time behavior of the 
quantum system in thermal medium allows the process such as x V-PX j even if 
it is kinematically forbidden for the on-shell S-matrix element. 

The essential part of our approach is that we separate the system (p part and 
the environment x part, and integrate over the environment part, assuming a ther- 
mal equilibrium for the lighter particle x- We thus implicitly assume that x has 
interaction with other light particles or among themselves not written in the above 
hamiltonian, to ensure the thermal equilibrium for x- 

Description of this class of the open system is conveniently done in the path 
integral approach, the influence functional method ||. The x integration is of a 
Gaussian type, including the thermal ensemble average, and one obtains a non-local 
exponent, symbolically of the form, (p 2 (l) a(l — 2) <^ 2 (2) in the influence functional. 
The kernel a is known and complex, reflecting the dissipation in the open system. 
We have developed a new Hartree type of approximation, replacing four ip's here by 
two ip's multiplied by the correlator G = i (<p<p); /3(1 , 2) = — i a(l — 2) G(l , 2) . The 
exponent of this Hartree model <p(l)/3(l , 2)<p(2) contains the two-body correlator 
G to be determined self-consistently. 

Derivation of the self-consistency equation is based on a physical assumption that 
time variation of the particle distribution function occurs more slowly than indivisual 
microscopic reactions occur. The two time argument in the correlator (<£>(1)<£>(2)) has 
then a slowing varying center of mass (CM) time (t± + t 2 )/2, which can be taken 
fixed for discussion of the short time dynamics. A similar separation of the CM and 
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the relative coordinate to distinguish the slow and the fast process is also taken in 
other approaches @, 0. 

As to the short time dynamics our reduced Hartree model is equivalent to a 
Gaussian model that consists of an infinite set of oscillators || defined by the given 
kernel above, once the CM time is fixed. The oscillator model in this form is exactly 
solvable and a convenient form is in ||. The time evolution, with regard to the 
relative time t\ —ti, exhibits an exponential decay law, and one may take the infi- 
nite relative time limit which has a non-vanishing equilibrium contribution for the 
relevant correlators. This equilibrium value is then allowed to slowly change with 
the CM time (t 1 + t 2 )/2. 

The coincident time limit of two-body correlators needed for our calculation of 
the occupation number 

,/A , <Pk<P-k U k 1 

/(«) = ( ~2=— + y WP-fc > - £ W 

(where Uk is a reference energy for the momentum k mode, taken to be the renor- 
malized u k (T), eq.(§), and <p k is the Fourier k mode of the p field) is calculable, 
for instance, in the operator approach H. The equilibrium occupation number thus 
calculated satisfies a self-consistency equation (CM time being omitted here); 

uk) + - = — r r+( ^ } _ — , ( 7) 

qV 2 4u k 7-oc (uj - u k {T)Y + (nr^uj, , k)/2uj k ) 2 
Ql(T) = P + Ml + A T 2 + 5Il(u k , k) . (8) 

Here r±(u , k) = r(u , k) ± r(— uo , k) , where (/3 = 1/T) 

( l\ o f d3k> ( r x (uj + uj k < ,k + k') - 

is the Fourier transform of the kernel f3(x\ , x 2 ) with respect to the relative coordinate 



x\ — %i- Here u k = J k 2 + M\. In this self-consistency equation for / eq the spectral 
function r x for the two-x state is a given function; it is odd in k and 

A 2 / 2 1 - e-^+l \ 

,k) = — e(k ) f 6(\k \ -k) + -ln 1 _ e _, | ,_ | J , (10) 
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with k = \k\ and u± = (\ko\ ± k)/2 , assuming that the x particle is massless. The 
temperature dependence of the ip mass is incorporated to 0[A 2 ], including the finite 
part of the proper self-energy SU after the mass (Mr being the renormalized ip mass, 
hereafter denoted by M for simplicity) and the wave function renormalization. 

The continuous energy integral flj]) has a Breit-Wigner form, at the peak of 
uj = u)fc(T) and of width T k = irr-(ui k , k)/u k . At high temperatures of T > M 
the entire Breit-Wigner peak region is integrated, giving the formula in the narrow 
width limit, f e q{k) ~ r(— u k ,k)/ r_(uj k ,k). At intermediate and low temperatures, 
a better approximation is given by a sum of the small resonance term plus a larger 
continuous integral; 

Uk) = ri ~ lJk '* ] +6Uk), (11) 



5/eq(fc) = -7— / ^ 7 F2/ . • 12 

4uj k J-oc (uj - uj k (T)y + T% A 



The tilded 5/ eq contains terms to be renormalized away by subtraction, the renor- 
malized one being denoted by 5f eq hereafter. 

A straightforward computation gives the off-shell contribution after renormaliza- 
tion; to the leading order of T/M 



6 fe q (k)^——— I dq 



167r 4 ku k Jo 2q + C(2)T e<?/ 2T - 1 \u; fc + u; fc _ ? w fc + u; fe+9/ 



(13) 



Physical processes that contribute to this result are predominantly inverse annihila- 
tion XX ~ * an d 1 to 3 process, x ~^ XV^V 9 ; which is only a small fraction of the 
entire contribution. The momentum integration of this occupation number gives 

c A 2 T 4 

Sn ea « — , (14) 

q 192tt 3 M K J 

oo t»2 i 

da; — ■ ^ 0.270. (15) 

The time evolution equation may be derived for the occupation number f(k). 
This equation is non-Markovian, containing an initial memory term. A simple 
Markovian approximation becomes possible, again under the assumption of the slow 
time variation of the <p number density. The Markovian equation thus derived takes 
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r fc (/(M)-/eq(M)) • (is) 



/ ^ { r t nn f(k,t)-rrUk,t)) , (17) 



the form, 

df(k,t) 
dt 

The usual Boltzmann equation follows when one approximates / eq by dropping 
the off-shell contribution 5f eq . One may use the explicit formula for (roughly 
annihilation cross section times target number density — inverse rate), and write the 
right side of it ( r_(uJk ,k)f(k) — r(— cu^ , fc) ) in the form equivalent to the Boltzmann 
equation in the thermal x medium. 

Integration of the distribution function over the momentum gives a rate equation 
for the number density which is of primary interest in the annihilation-scattering 
problem. By the momentum integration the scattering contribution nearly drops 
out. This is reasonable, because the scattering does conserve the ip particle number, 
hence the scattering process does not cause the change of the p number density. 

The evolution equation for the number density is simplified considering the can- 
cellation of scattering terms, to give 

dn r d 3 k 
~dt ~ ~~ J (2?r) 

where r^ nn (T™ v ) is the rate keeping the annihilation (inverse annihilation) term. 
The first destructive term in the right side gives the usual, thermally averaged 
rate (cr a v) n 2 , while the second production term gives, at low temperatures, ~ 
{<?av) n th 5n eq , with n th = ((3)T 3 /n 2 and 5n eq given by (|TJ). 

In cosmology the thermal environment gradually changes according to the adia- 
batic law; the temperature cx the cosmic scale factor l/a(t). Along with the number 
density oc l/a 3 (t), the time derivative operator is modified to % + 3 if n = ■ ■ ■ . 
Here H = a/a is the Hubble parameter. 

A detailed behavior of the ip number density may be worked out by examining 

^ + 3Hn = - M ( n 2 - 6^ T s - < B ) , (18) 

where »mb ~ (MT/2n) 3 / 2 e~ M l T is the thermal number density of zero chemical 
potential. The evolution equation is further simplifed by = T~ 3 ( % + 3H n ) . 
The equation thus derived is our main result stated in the first part of this paper for 
the S-wave annihilation. The freeze-out parameters in the boson-pair annihilation 
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model become T f /M « 700 (N^/X 2 ) (M/m pl ) 2 / 3 , Y f « 0.06 JV 1 ^ (M/m^Va , 
thus giving the closure density atM«l GeV. 
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Figure caption 

Fig.l 

Time evolution for the number density divided by temperature 3 . Parameters 
are taken from the boson-pair annihilation model in the text; n = and 5/i] = 
4.6 x 1(T 4 VN M/rripx , 5 = 5.5 x KT 6 A 2 . The case for M = 1 GeV , N = 43/4 and 
indicated values of A is shown, along with the evolution when the on-shell Boltzmann 
equation is used. 

Fig.2 

Contour map for the present mass density equated to the closure p c = 1 x 
10 4 eV cm~ 3 , in the parameter space (S/rj , assuming the relation of the boson 
model, S/rj = 1.2 x 10 -22 M/GeV . The contours for 0.1 p c and 0.01p c are also 
shown. 
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Fig 1 
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